/********************************************************************************
*  This file is part of CinoLib                                                 *
*  Copyright(C) 2016: Marco Livesu                                              *
*                                                                               *
*  The MIT License                                                              *
*                                                                               *
*  Permission is hereby granted, free of charge, to any person obtaining a      *
*  copy of this software and associated documentation files (the "Software"),   *
*  to deal in the Software without restriction, including without limitation    *
*  the rights to use, copy, modify, merge, publish, distribute, sublicense,     *
*  and/or sell copies of the Software, and to permit persons to whom the        *
*  Software is furnished to do so, subject to the following conditions:         *
*                                                                               *
*  The above copyright notice and this permission notice shall be included in   *
*  all copies or substantial portions of the Software.                          *
*                                                                               *
*  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR   *
*  IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,     *
*  FITNESS FOR A PARTICULAR PURPOSE AND NON INFRINGEMENT. IN NO EVENT SHALL THE *
*  AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       *
*  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      *
*  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS *
*  IN THE SOFTWARE.                                                             *
*                                                                               *
*  Author(s):                                                                   *
*                                                                               *
*     Marco Livesu (marco.livesu@gmail.com)                                     *
*     http://pers.ge.imati.cnr.it/livesu/                                       *
*                                                                               *
*     Italian National Research Council (CNR)                                   *
*     Institute for Applied Mathematics and Information Technologies (IMATI)    *
*     Via de Marini, 6                                                          *
*     16149 Genoa,                                                              *
*     Italy                                                                     *
*********************************************************************************/
#include <cinolib/laplacian.h>
#include <cinolib/symbols.h>
#include <Eigen/Sparse>

namespace cinolib
{

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

template<class M, class V, class E, class P>
CINO_INLINE
std::vector<Eigen::Triplet<double>> laplacian_matrix_entries(const AbstractMesh<M,V,E,P> & m,
                                                             const int mode,
                                                             const int n) // diagonally replicate n times
{
    std::vector<Entry> entries;

    uint nv = m.num_verts();
    uint base[n];
    for(int i=0; i<n; ++i) base[i] = nv*i;

    for(uint vid=0; vid<m.num_verts(); ++vid)
    {
        std::vector<std::pair<uint,double>> wgts;
        m.vert_weights(vid, mode, wgts);
        double sum = 0.0;
        for(auto item : wgts)
        {
            for(int i=0; i<n; ++i)
            {
                entries.push_back(Entry(base[i] + vid, base[i] + item.first, item.second));
            }
            sum -= item.second;
        }
        if(sum == 0.0)
        {
            std::cerr << "WARNING: null row in the matrix! (disconnected vertex? I put 1 in the diagonal)" << std::endl;
            sum = 1.0;
        }
        for(int i=0; i<n; ++i)
        {
            entries.push_back(Entry(base[i] + vid, base[i] + vid, sum));
        }
    }

    return entries;
}

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

template<class M, class V, class E, class P>
CINO_INLINE
Eigen::SparseMatrix<double> laplacian(const AbstractMesh<M,V,E,P> & m, const int mode, const int n)
{
    std::vector<Entry> entries = laplacian_matrix_entries(m, mode, n);

    uint nv = n*m.num_verts();
    Eigen::SparseMatrix<double> L(nv,nv);
    L.setFromTriplets(entries.begin(), entries.end());

    return L;
}

}
